###################################################################################################
## Payson (2020) "The Partisan Logic of City Mobilization" Replication Code
## Panel Data Analysis (Table 1 and Figures 2 - 3)
## Created with R version 3.5.3 (2019-03-11) -- "Great Truth"
###################################################################################################

## Set working directory
setwd()


## Load pacakges
library(ggplot2); library(reshape); library(effects); library(lmtest); 
library(sandwich); library(stargazer); library(lfe); library(dplyr)


## Functions for tables (Note: tables will not compile without these functions)
pretty <- function(x) {
  paste(prettyNum(x, big.mark = ","))}

num.cities <- function(x) {
  length(which(getfe(x)=="fips"))}

N <- function(x) {
  x$N}

getmean <- function(x) {
  round(mean(x$response), 2)}


###################################################################################################

## load panel.csv
dta <- read.csv("panel.csv")



###################################################################################################
## Table 1: Effect of Partisan Mismatch on City Lobbying
###################################################################################################

## There were 732 elections resulting in a partisan mismatch between cities and their state reps
table(dta$mismatch)

m <- felm(lobby ~ mismatch + log.pop | fips + year | 0 | fips, 
          data = subset(dta, pres.dem.tercile != "Moderate"))

m1 <- felm(lobby ~ mismatch + log.pop | fips + state.by.year | 0 | fips, 
           data = subset(dta, pres.dem.tercile != "Moderate"))


m2 <- felm(lobby ~ mismatch + log.pop + log.income + log.own.source +
             pct.white + log.med.house | fips + state.by.year | 0 | fips, 
           data = subset(dta, pres.dem.tercile != "Moderate"))


## Table 1: Effect of Partisan Mismatch on City Lobbying
stargazer(m, m1, m2, omit.table.layout = "s", type = "text",
          star.cutoffs = c(0.05), keep = c("mismatch", "log.pop"),
          title = c("Effect of Partisan Mismatch on City Lobbying."),
          dep.var.caption = "Probability of Lobbying", dep.var.labels.include = FALSE,
          covariate.labels = c("Partisan Mismatch", "Population (Log)"),
          add.lines = list(c("City FEs", "Y", "Y", "Y"),
                           c("Year FEs", "Y", ""),
                           c("State-Year FEs", "", "Y", "Y"),
                           c("Full Controls", "", "", "Y"),
                           c("Observations", pretty(m$N), pretty(m1$N), pretty(m2$N)),
                           c("\\# Cities", pretty(num.cities(m)), 
                             pretty(num.cities(m1)), pretty(num.cities(m2))),
                           c("Mean Lobbying Probability", getmean(m), getmean(m1), 
                             getmean(m2))),
          notes.append = FALSE, notes.align = "l", notes.label  = "",
          notes = "Robust standard errors clustered by city. $^{*}$p$<$0.05")



##################################################################################################
## Figure 1: Partisan Mismatch: Leads and Lags
###################################################################################################

m <- felm(lobby ~  mismatch.lead2 + mismatch.lead1 + mismatch + mismatch.lag1 + log.pop + 
            log.income + log.own.source + pct.white | fips + state.by.year | 0 | fips, 
          data = subset(dta, pres.dem.tercile != "Moderate"))

## Format data for figure
beta <- c(m$coefficients[1], m$coefficients[2], m$coefficients[3])
se <- c(m$cse[1], m$cse[2], m$cse[3])
figdata <- data.frame(beta = c(beta), se = c(se), x = c(-.5, 0, .5))

## Figure: Leads and Lags
(p <- ggplot(figdata, aes(x = x, y = beta)) + 
    geom_point(size = 2) +
    geom_errorbar(aes(ymin = beta - (1.9 * se), ymax = beta + (1.9 * se)), 
                  size = .5, width = .1) +
    geom_hline(yintercept = 0, linetype = 2) +
    scale_x_continuous(breaks = c(-.5, 0, .5), 
                       labels = c("2 Years Prior", "1 Year Prior", "Mismatch")) +
    geom_vline(xintercept = .25) +
    theme_bw(base_family = "serif", base_size = 12)  +
    labs(y = "Effect on Lobbying Probability", 
         x = "Years Relative to Treatment"))



###################################################################################################
## Figure 2: Correlation Between City and State Representative Ideology
###################################################################################################

year2007 <- subset(dta, year == 2007)

NY <- year2007[year2007$city == "NEW YORK CITY", ]
Berkeley <- year2007[year2007$city == "BERKELEY CITY" & year2007$state=="California", ]
Amarillo <- year2007[year2007$city == "AMARILLO CITY", ]


## Figure: Correlation Between City and State Representative Ideology
(p <- ggplot(year2007, aes(x = mrp.estimate, y = mean.np.lower)) +
  geom_point(alpha = .4) +
  xlim(-1.1, .7) +
  stat_smooth(method = "lm", colour = "black") +
  labs(y = "State Representative Ideology", x = "City Ideology") +
  theme_bw(base_family = "serif", base_size = 13) +
  theme(plot.background = element_blank(),
        panel.grid.major = element_blank(),
        panel.grid.minor = element_blank()) +
  annotate("text", x = NY$mrp.estimate, y = NY$mean.np.lower, 
           label = "New York City", size = 4.5) +
  annotate("text", x = Berkeley$mrp.estimate + .1, y = Berkeley$mean.np.lower - .1, 
           label = "Berkeley, CA", size = 4.5) +
  annotate("text", x = Amarillo$mrp.estimate + .05, y = Amarillo$mean.np.lower - .4, 
           label = "Amarillo, TX", size = 4.5))



###################################################################################################
##Figure 3: Effect of Ideological Distance on Lobbying
###################################################################################################

## Note: be sure to set correct levels on np.quantile variable prior to analysis
dta$np.quantile <- factor(dta$np.quantile, levels = c("Most Liberal", "Liberal", "Moderate", 
                                                       "Conservative", "Most Conservative"))

## Figure: Effect of Ideological Mismatch on Lobbying
m <- lm(lobby.lead1 ~ np.quantile + log.pop + log.own.source + log.income +
          pct.white + log.med.house + fips + state.by.year, 
        data = subset(dta, mrp.tercile == "Liberal"))


m1 <- lm(lobby.lead1 ~ np.quantile + log.pop + log.own.source + log.income +
           pct.white + log.med.house + fips + state.by.year, 
         data = subset(dta, mrp.tercile == "Conservative"))


## Use ?effects for additional information about the effects pacakges used to
## create the effects plots
plot.dta <- data.frame(rep.type = rep(c(0, 2, 4, 6, 8)),
                       city.type = c(rep("Liberal City", 5), rep("Conservative City", 5)),
                       lower = c(summary(effect("np.quantile", m))$lower[1],
                                 summary(effect("np.quantile", m))$lower[2],
                                 summary(effect("np.quantile", m))$lower[3],
                                 summary(effect("np.quantile", m))$lower[4],
                                 summary(effect("np.quantile", m))$lower[5],
                                 summary(effect("np.quantile", m1))$lower[1],
                                 summary(effect("np.quantile", m1))$lower[2],
                                 summary(effect("np.quantile", m1))$lower[3],
                                 summary(effect("np.quantile", m1))$lower[4],
                                 summary(effect("np.quantile", m1))$lower[5]),
                       upper = c(summary(effect("np.quantile", m))$upper[1],
                                 summary(effect("np.quantile", m))$upper[2],
                                 summary(effect("np.quantile", m))$upper[3],
                                 summary(effect("np.quantile", m))$upper[4],
                                 summary(effect("np.quantile", m))$upper[5],
                                 summary(effect("np.quantile", m1))$upper[1],
                                 summary(effect("np.quantile", m1))$upper[2],
                                 summary(effect("np.quantile", m1))$upper[3],
                                 summary(effect("np.quantile", m1))$upper[4],
                                 summary(effect("np.quantile", m1))$upper[5]),
                       effect = c(summary(effect("np.quantile", m))$effect[1],
                                  summary(effect("np.quantile", m))$effect[2],
                                  summary(effect("np.quantile", m))$effect[3],
                                  summary(effect("np.quantile", m))$effect[4],
                                  summary(effect("np.quantile", m))$effect[5],
                                  summary(effect("np.quantile", m1))$effect[1],
                                  summary(effect("np.quantile", m1))$effect[2],
                                  summary(effect("np.quantile", m1))$effect[3],
                                  summary(effect("np.quantile", m1))$effect[4],
                                  summary(effect("np.quantile", m1))$effect[5]))


limits <- aes(ymin = lower, ymax = upper)

(p <- ggplot(plot.dta, aes(rep.type, effect)) +
    geom_point() +
    facet_wrap(~city.type) +
    geom_line() + ylim(.2, .8) +
    scale_x_continuous(labels = c("Most \n Liberal", "Liberal", "Moderate", 
                                  "Conservative", "Most \n Conservative \n")) +
    geom_errorbar(limits, width = 0.1) +
    theme_bw(base_family = "serif", base_size = 16) +
    theme(panel.spacing = unit(2, "lines"),
          plot.margin = unit(c(1,2,1,1), "cm"),
          plot.background = element_blank(),
          strip.background = element_rect(fill = "white")) +
    labs(y = "Probability of Lobbying", x = "State Representative Ideology"))
